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Abstract. The enhancement and detection of elongated structures in 
noisy image data is relevant for many biomedical applications. To han¬ 
dle complex crossing structures in 2D images, 2D orientation scores 
U : R 2 x S 1 —y R were introduced, which already showed their use 
in a variety of applications. Here we extend this work to 3D orientation 
scores U : R 3 x S' 2 —>> R. First, we construct the orientation score from a 
given dataset, which is achieved by an invertible coherent state type of 
transform. For this transformation we introduce 3D versions of the 2D 
cake-wavelets, which are complex wavelets that can simultaneously de¬ 
tect oriented structures and oriented edges. For efficient implementation 
of the different steps in the wavelet creation we use a spherical harmonic 
transform. Finally, we show some first results of practical applications of 
3D orientation scores. 

Key words: Orientation Scores, Reproducing Kernel Spaces, 3D Wavelet 
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1 Introduction 

The enhancement and detection of elongated structures is important in many 
biomedical image analysis applications. These tasks become problematic when 
multiple elongated structures cross or touch each other in the data. In these 
cases it is useful to decompose an image in local orientations by constructing an 
orientation score. In the orientation score, we extend the domain of the data to 
include orientation in order to separate the crossing or touching structures (Fig. 
1). From 3D data / : M 3 —>> M we construct a 3D orientation score U : M 3 x S 2 —» 
M, in a similar way as is done for the more common case of 2D data / : M 2 —>• M 
and 2D orientation score U : M 2 x S 1 M. Next, we consider operations on 
orientation scores, and process our data via orientation scores (Fig. 2). For such 
operations it is important that the orientation score transform is invertible, in a 
well-posed manner. In comparison to continuous wavelet transforms on the group 
of 3D rotations, translations and scalings, we use all scales simultaneously and 
exclude the scaling group from the wavelet transform and its adjoint, yielding a 
coherent state type of transform [1], see App.A. This makes it harder to design 
appropriate wavelets, but has the computational advantage of only needing a 
single scale transformation. 
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The 2D orientation scores have already showed their use in a variety of 
applications. In [11,17] the orientation scores were used to perform crossing¬ 
preserving coherence-enhancing diffusions. These diffusions greatly reduce the 
noise in the data, while preserving the elongated crossing structures. Next to 
these generic enhancement techniques, the orientation scores also showed their 
use in retinal vessel segmentation [3], where they were used to better handle 
crossing vessels in the segmentation procedure. 

To perform detection and enhancement operations on the orientation score, 
we first need to transform a given greyscale image or 3D dataset to an orientation 
score in an invertible way. In previous works various wavelets were introduced to 
perform a 2D orientation score transform. Some of these wavelets did not allow 
for an invertible transformation (e.g. Gabor wavelets [15]). A wavelet that allows 
an invertible transformation was proposed by Kalitzin [14]. A generalization 
of these wavelets was found by Duits [8] who derived a unitarity result and 
expressed the wavelets in a basis of eigenfunctions of the harmonic oscillator. 
This type of wavelet was also extended to 3D. This wavelet however has some 
unwanted properties such as poor spatial localization (oscillations) and the fact 
that the maximum of the wavelet did not lie at its center [8, Fig. 4.11]. In [8] 
a class of cake-wavelets were introduced, that have a cake-piece shaped form in 
the Fourier domain (Fig. 5). The cake-wavelets simultaneously detect oriented 
structures and oriented edges by constructing a complex orientation score U : 
M 2 x S 1 —> C. Because the different cake-wavelets cover the full Fourier spectrum, 
invertibility is guaranteed. 

In this paper we propose an extension of the 2D cake-wavelets to 3D. First, 
we discuss the theory of invertible orientation score transforms. Then we con¬ 
struct 3D cake-wavelets and give an efficient implementation using a spherical 
harmonic transform. Finally we mention two application areas for 3D orienta¬ 
tion scores and show some preliminary results for both of them. In the first 
application, we present a practical proof of concept of a natural extension of 
the crossing preserving coherence enhancing diffusion on invertible orientation 
scores (CEDOS) [11] to the 3D setting. Compared to the original idea of co¬ 
herence enhancing diffusion acting directly on image-data [18,4,5] we have the 
advantage of preserving crossings. Diffusions on SE(3) have been studied in pre¬ 
vious SSVM-articles, see e.g. [6], but the full generalization of CEDOS to 3D 
was never established. 

2 Invertible Orientation Scores 

An invertible orientation score [/] : M 3 x S 2 C is constructed from a given 
ball-limited 3D dataset / G L^M 3 ) = {/ G L 2 (M 3 )|supp(J 7 /) C with 

g > 0 by correlation * with an anisotropic kernel 

(Wb[/D( X ^ n ) = C0n*/)(x) = f ^n(x' -x)/(x')dx / , (1) 

J R 3 

where G L 2 (M 3 ) DLi(M 3 ) is a wavelet aligned with and rotationally symmetric 
around the 2 -axis, and ^ n (x) = ^(R^x) G L 2 (M 3 ) the rotated wavelet aligned 
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with n. Here R n is any rotation which rotates the 2 -axis onto n where the specific 
choice of rotation does not matter because of the rotational symmetry of 'ip. The 
over line denotes a complex conjugate. The exact reconstruction formula for this 
transformation is 


/(x) = (wy[w4/]])(x) 



x 1 —y 


[ ( 

Js 2 


C0n* Wv, [/](■, n))(x) dcr(n) 


( x ), 


( 2 ) 


with J ^3 the Fourier transform on M 3 given by (F f)( uj) = (2i r) 2 J R3 e 2 u; ' x /(x)dx 
and ^ n (x) = '0 n (—x). In fact is a unitary mapping on to a reproducing ker¬ 
nel space, see App. A. The function M ^ : M 3 —)> M + is given by 


M^(u) = (27t) 2 f |J’ R3 ['0n](^)| 2 dcr(n). (3) 

Js 2 

The function M ^ quantifies the stability of the inverse transformation [8], since 
specifies how well frequency component u; is preserved by the cascade 
of construction and reconstruction when M^ 1 would not be included in Eq. (2). 
An exact reconstruction is possible as long as 


3m>o,s>o 0 < S < < M < 00 , for all u = B$, e . (4) 

In practice it is best to aim for « 1, in view of the condition number of 

: L^M 3 ) L^M 3 x S 2 ) with W^f = Also, when = 1 we have 

L 2 -norm preservation 

II/IIl 2 (R3) = II^/IIl 2 (R3 xS2) , for all / € Lf(M 3 ), (5) 

and Eq. ( 2 ) simplifies to /(x) = / s2 (^ n * W- 0 [/](•, n))(x)dcr(n). We can further 
simplify the reconstruction for wavelets for which (27r)^ f s2 [xp n ](u;)da(n) ~ 
1 , where the reconstruction formula simplifies to an integration over orientations 


/(x) ~ [ VV^/(x, n) dcr(n). 

Js 2 


( 6 ) 




Fig. 1. 2D Orientation score for an exemplary im- Fig. 2. A schematic view of 
age. In the orientation score crossing structures are image processing via invertible 
disentangled because the different structures have a orientation scores, 
different orientation. 
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/(x) V>e x (x) HV[/](x,e x ) 



V»n, (x) Wy, [/] (x, m) Vfy [/] (x, n) 



Fig. 3. Creating a 3D orientation score. Top: The data / is correlated with an oriented 
filter f i/j ex to detect structures aligned with the filter orientation e x . Bottom left: This 
is repeated for a discrete set of filters with different orientations. Bottom right: The 
collection of 3D datasets constructed by correlation with the different filters is an 
orientation score and is visualized by placing a 3D dataset on a number of orientations. 


2.1 Discrete Invertible Orientation Score Transformation 


In the previous section, we considered a continuous orientation score transfor¬ 
mation. In practice, we have only a finite number of orientations. To determine 
this discrete set of orientations we uniformly sample the sphere using platonic 
solids and/or refine this using tessellations of the platonic solids. 

Assume we have a number N 0 of orientations V = {ni, ri 2 ,..., rijv 0 } C -S' 2 , 
and define the discrete invertible orientation score W^[f] : M 3 x V C by 


(W$[/])(x,ni) = On ; */)(x). (7) 

The exact reconstruction formula is in the discrete setting given by 
/(x) = ((W^) -1 [W^ [/]]) (x) 




— 1 
R 3 




N 0 


x-> 5^(^n i *W^[/](-,n i ))(x)dcr(n i ) 


i =1 


( x )> 


(8) 
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with da(rii) the discrete spherical area measure which for reasonably uniform 
spherical sampling can be approximated by dcr(n^ ) « ft, and 

No 

= (2tt) 5 y]|J’ R 3[^ n J(w)| 2 d<j(n i ). (9) 

i= 1 

Again, an exact reconstruction is possible iff 0 < S < M$(u) < M < oo. 

3 3D Cake-Wavelets 

A class of 2 D cake-wavelets, see [ 8 ], was successfully used for the 2 D orientation 
score transformation. We now generalize these 2D cake-wavelets to 3D cake- 
wavelets. Our 3D transformation using the 3D cake-wavelets should fulfill a set 
of requirements, compare [ 11 ] : 

1 . The orientation score should be constructed for a finite number (N 0 ) of 
orientations. 

2. The transformation should be invertible and all frequencies should be trans¬ 
ferred equally to the orientation score domain ( M^ « 1 ). 

3. The kernel should be strongly directional. 

4. The kernel should be polar separable in the Fourier domain, i.e., ( J = 
g(p)h{0 ,0), with uj = (uj x ,uj y ,uj z ) = (p sin # cos 0, p sin # sin </>, pcos#). Be¬ 
cause by definition the wavelet ^ has rotational symmetry around the z-axis 
we have h(0, <p) = 

5. The kernel should be localized in the spatial domain, since we want to pick 
up local oriented structures. 

6 . The real part of the kernel should detect oriented structures and the imag¬ 
inary part should detect oriented edges. The constructed oriented score is 
therefore a complex orientation score. 


3.1 Construction of Line and Edge Detectors 

We now discuss the procedure used to make 3D cake-wavelets. According to 
requirement 4 we only consider polar separable wavelets in the Fourier domain, 
so that tc?) = p(p)fj(0). For the radial function g(p) we use, as in [11], 

2 N ( n 2 f -l\k 

g(p) = M N (p 2 t~ l ) = e~^r ^ , ( 10 ) 

k =0 

which is a Gaussian function with scale parameter t multiplied by the Taylor 
approximation of its reciprocal to order N to ensure a slower decay. This func¬ 
tion should go to 0 when p tends to the Nyquist frequency p^. Therefore the 
inflection point of this function is fixed at ypw with 0 < 7 < 1 by setting 
t = 2 ^+2N • P ra( Aice we have g = p^v, and because radial function g causes 
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to become really small when coming close to the Nyquist frequency, recon¬ 
struction Eq.(8) becomes unstable. We solve this by either using approximate 
reconstruction Eq.(6) or by replacing M^ —)• max(M^,e), with e small. Both 
make the reconstruction stable at the cost of not completely reconstructing the 
highest frequencies which causes some additional blurring. 

We now need to find an appropriate angular part f) for the cake-wavelets. 
First, we specify an orientation distribution A : S 2 —>• M + , which determines 
what orientations the wavelet should measure. To satisfy requirement 3 this 
function should be a localized spherical window, for which we propose a B-spline 
A(6, (f)) = with sq > 0 and B k the fcth order B-spline given by 


B k {x) = (B k ~ x * B°)(x), B°{x) 


1 if — \ < x < \ 
0 otherwise 


( 11 ) 


The parameter sq determines the trade-off between requirements 2 and 3, where 
higher values give a more uniform M^ at the cost of less directionality. 

First consider setting h = A so that ^ has compact support within a convex 
cone in the Fourier domain. The real part of the corresponding wavelet would 
however be a plate detector and not a line detector (Fig. 4). The imaginary part 
is already an oriented edge detector, and so we set 


= A(0,4>) -A(ir-O,0 + ir) = B k (T) - B k , (12) 

where the real part of the earlier found wavelet vanishes by anti-symmetrization 
of the orientation distribution A while the imaginary part remains. As to the 
construction of hR e , there is the general observation that we detect a structure 
that is perpendicular to the shape in the Fourier domain, so for line detection 
we should aim for a plane detector in the Fourier domain. To achieve this we 
apply the Funk transform to A , and we define 

h Re (d, 0) = FA(9, </>)=[ A(n') ds(n'), (13) 

JS p (n(0,<l>)) 

where integration is performed over S p ( n) denoting the great circle perpendic¬ 
ular to n(0, 4b) = (sin 0 cos 0, sin 6 sin 0, cos 0). This transformation preserves the 
symmetry of A , so we have hR e (0 , 4>) = §Re(Q)- Thus, we finally set 

m = i)Re(0) + I Mm{0). (14) 

For an overview of the transformations see Fig. 5. 

3.2 Efficient Implementations Via Spherical Harmonics 

In Subsection 3.1 we defined the real part and the imaginary part of the wavelets 
in terms of a given orientation distribution. In order to efficiently implement the 
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various transformations (e.g. Funk transform), and to create the various rotated 
versions of the wavelet we express our orientation distribution A in a spherical 
harmonic basis {Yj m } up to order L: 

l i 

A{9, <£) = £ £ c ^mYr(0, </>), Le N. (15) 

1=0 m=—l 

Because of the rotational symmetry around the 2 )-axis, we only need the spher¬ 
ical harmonics with m = 0, i.e., A(6,(ft) = Q,oFJ°($, (ft)- For determining 

the spherical harmonic coefficients we use the pseudo-inverse of the discretized 
inverse spherical harmonic transform (see [9, Section 7.1]), with discrete orien¬ 
tations given by an icosahedron of tesselation order 15. 

Funk Transform According to [7], the Funk transform of a spherical harmonic 
equals 

FYr(e,4>)= f Y l m (n')ds(n') = 2irP l (0)Y l m (e,(f>), (16) 

J S p (n(d,(f>)) 

with P/(0) the Legendre polynomial of degree l evaluated at 0. We can therefore 
apply the Funk transform to a function expressed in a spherical harmonic basis 
by a simple transformation of the coefficients c™ — 27rPi(0)c] n . 

Anti-Symmetrization We have Y l m (7r — 0, (ft A tt) = (—l) z YJ m (0, (ft). We there¬ 
fore anti-symmetrize the orientation distribution Eq. (12) via c f 1 —>• (1 — (—l) z )c[ n . 

Making Rotated Wavelets To make the rotated versions of wavelet ^ we 
have to find h n in = g(p)h n (6 , (ft). To achieve this we use the steerability of the 
spherical harmonic basis. Spherical harmonics rotate according to the irreducible 
representations of the SO(3) group D l m m ,(a,/3, 7 ) (Wigner-D functions) 

1 

Kn^,Yr(0,4>) = £ D l m!m ,(a,p, 7 )Yr\e^). (17) 

m'=l 

Here a,/3 and 7 denote the Euler angles with counterclockwise rotations, i.e., 
R — Re z ,/3-f^e z ,7 • This gives 

L l l 

MM) = = £ £ £ a^Dl^p^Yr'ie^). (18) 

1=0 m=—lm' =—l 

Because both anti-symmetrization and Funk transform preserve the rotational 
symmetry of A , we have h(0, (ft) = Yl^=o &z,oFJ°(0, (ft), and Eq. (18) reduces to 

l 1 

h n (9,<j)) = £ £ a lfi D l Otm ,(0,^)Yr'(e,cP). 

1=0 mJ =—l 


(19) 
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A(0,<t>) (Tip){w)=g(p)A(e,<j>) Im[ip\(x) ife#](x) 



Fig. 4. When directly setting orientation distribution A as angular part of the wavelet 
h we construct plate detectors. From left to right: Orientation distribution A, wavelet 
in the Fourier domain, the plate detector (real part) and the edge detector (imaginary 
part). Orange: Positive iso-contour. Blue: Negative iso-contour. Parameters used: L = 
16, sq = 0.6, k = 2,N = 20,7 — 0-8 and evaluated on a grid of 51x51x51 pixels. 



Fig. 5. Cake-Wavelets. Top: 2D cake-wavelets. From left to right: Illustration of the 
Fourier domain coverage, the wavelet in the Fourier domain and the real and imag¬ 
inary part of the wavelet in the spatial domain. [3]. Bottom: 3D cake-wavelets. 
Overview of the transformations used to construct the wavelets from a given ori¬ 
entation distribution. Upper part: The wavelet according to Eq. (12). Lower part: 
The wavelet according to Eq. (13). IFT: Inverse Fourier Transform. Parameters used: 
L = 16, sq = 1.05, k = 2,N = 20,7 — 0-8 and evaluated on a grid of 31x31x31 pixels. 






























Invertible Orientation Scores of 3D Images 9 


4 Applications 

4.1 Adaptive Crossing Preserving Flows 

We now use the invertible orientation score transformation to perform data- 
enhancement according to Fig. 2 . Because M 3 x s 2 is not a Lie group, it is 
common practice to embed the space of positions and orientations in the Lie 
group of positions and rotations SE(3) by setting 

U (x, R) = E/(x, R • e*), U (x, n) = U (x, R n ), (20) 

with R n any rotation for which R n • e z = n. This holds in particular for orien¬ 
tation scores U = Wp/. The operations <3? which we consider are scale spaces on 
SE(3) (diffusions), and are given by = <f>t with 

MU)(y,n)=W(y,R n ,t). ( 21 ) 

Here W is the solution of 

dW ® - -—-— 

— (5,0 = £ AilgDijAjlgWM, W\t=o = WM (22) 
i,j =1 

where in coherence enhancing diffusion on orientation scores (CEDOS) Dij is 

adapted locally to data W^[f] based on exponential curve fits (see [10]), and 
with Ai| p= ( x> R) = (L g )*Ai | e the left-invariant vector fields on SE(3), for mo¬ 
tivation and details see [9]. Furthermore is chosen such that equivalence 
relation Eq. (20) is maintained for W. These operations are already used with¬ 
out adaptivity in the field of diffusion weighted MRI, where similar data (of the 
type M 3 x S 2 —>> M + ) is enhanced [9]. We then obtain Euclidean invariant image 
processing via 


Tf = w*/ xt O 0 O w+f = w; O P '40 O Wtf (23) 

which includes inherent projection P^, of orientation scores, even if @ = <& t maps 
outside of the space of orientation scores in the embedding space (see App. 
A). Below we show some preliminary results of these flows that enhance the 
elongated structures while preserving the crossing, Fig. 6 and Fig. 7. 


4.2 3D Vessel Tracking in Magnetic Resonance Angiography 
(MRA) Data 

We use the 3D orientation scores to extend the earlier work on 2D vessel segmen¬ 
tation via invertible orientation scores [3] to 3D vessel segmentation in MRA- 
data. Even though true crossing structures hardly appear in 3D data, we do 
encounter vessels touching other vessels/structures. The orientation scores also 
allow us to better handle complex structures, such as bifurcations. In Fig. 8 we 
show some first results of the vessel segmentation algorithm. 
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Fig. 6. Adaptive Crossing Preserving Flows. From left to right 3D visualization of 
artificial data, slice of data, slice of (data + Gaussian noise), slice of enhanced data. 
For the orientation score transformation we use: No = 42, so = 0.7, k = 2, N = 20, 7 = 
0.85, L = 16 evaluated on a grid of 21 x 21 x 21 pixels. We use approximate reconstruction 
Eq.( 8 ), and for diffusion we set t = 10. For the choice of Dij in CEDOS, see [10]. 



Fig. 7. Adaptive Crossing Preserving Flows combined with soft thresholding 
<£(£/) (x, n) = |C7(x, n)| 1-5 sgn(C7(x, n)) on data containing the Adam Kiewitzc vessel. 
From left to right: Slice of data, data after soft thresholding, data after CEDOS, data 
after CEDOS followed by soft thresholding. For parameters see Fig. 6 , but now t — 5. 




Fig. 8 . MRA vessel segmentation via in¬ 
vertible orientation scores. 


5 Conclusion 

We have extended 2D cake-wavelets to 3D cake-wavelets, which can be used 
for a 3D invertible orientation score transformation. Efficient implementation 
for calculating the wavelets via spherical harmonics were introduced. The devel¬ 
oped transformation allows us to consider all kinds of enhancement operations 
via orientation scores such as the adaptive crossing preserving flows which we 
are currently working on. Next to data-enhancement we also showed some first 
results of 3D vessel segmentation using 3D orientation scores. 
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A Invertible Orientation Scores of 3D-images and 
Continuous Wavelet Theory 


The continuous wavelet transform constructed by unitary irreducible represen¬ 
tations of locally compact groups was first formulated by Grossman et al. [13]. 
Given a Hilbert space H and a unitary irreducible representation g U g of any 
locally compact group G in H, a non-zero vector ^ £ H is called admissible if 


C* 


L 


[ 2 


d^ G (g) < oo, 


(24) 


where hq denotes the left-invariant Haar measure. Given an admissible vector 
^ and a unitary representation of a locally compact group G in H, the Coherent 
State (CS) transform : H —L 2 (G) is given by (W^[f])(g) = (U g ip, /)#. W ^ 
is an isometric transform onto a unique closed reproducing kernel space C^- 
with K^(g. g') = -^■{U g 'ip,U g f'ip)H as an L 2 -subspace [1]. 

We distinguish between the isometric wavelet transform : L^M 3 ) 
L 2 (G) and the unitary wavelet transform : L^M 3 ) —>• C^-. We drop the 
formal requirement of U being square-integrable and pj being admissible in the 
sense of (24), and replace the requirement by (4), as it is not strictly needed in 
many cases. This includes our case of interest G = SE( 3) and its left-regular 
action on L 2 (M 3 ) where gives rise to an orientation score : M 3 x S 2 —> C 


W^f(x,n) = W^f(x,R ri ), 


(25) 


with R n any rotation mapping e z onto n and ^ symmetric around the 2 -axis. 
Here the domain is the coupled space of positions and orientations: M 3 xi S 2 := 
S£7(3)/({0}xSO(2)),cf.[9]. 

From the general theory of reproducing kernel spaces [8, Thm 18],[2] (where 
one does not even rely on the group structure), it follows that : L^M 3 ) —)> 

US’S \y\ c<2 ]rp3 vi Cf2 

C K is unitary, where C K denotes the abstract complex reproducing ker¬ 
nel space consisting of functions on M 3 x S 2 with reproducing kernel 

■^f(y,n) (y 5 n ) = (^(y,R n )'0:^(y / ,R n /) V ; )l 2 (M 3 ) 5 (26) 

with left-regular representation (y, R) i—> given by (W( y? r/jVOM = 

^(R T (x — y)). Now, as the characterization of the inner product on 
is awkward [16], we provide a basic characterization next via the so-called M ^ 
inner product. This is in line with the admissibility conditions in [12]. 

Theorem 1. Let ip be such that (4) holds. Then : L^M 3 ) —>> is 

unitary, and we have 


(/^)l 2( m3) = (W^/,W^)m^ (27) 

where = (7m^ [>%,/], Tm^ with [Tm^[U}\{ y, n) : 


J- 


OJ i —¥ 


(27r)- 3 / 4 M- 1/2 ( W )^[C/(-,n)]( W ) 


(y)- 
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Proof. We rely on [17, Thm 1], where we set H = L 2 (R 3 ). The rest follows by 
well posed restriction to the quotient M 3 x s 2 . 

Corollary 1. Let M/ > 0 on M 3 . The space is a closed subspace of 

Hilbert space ® L 2 (5' 2 ), where = {/ G L 2 (M 3 )| M ^ 2 H[f] G L 2 (M 3 )}, 
and projection of embedding space onto the space of orientation scores is given 
by (P^(E/))(y,n) = (K (n)y) , t/) M * = (W^W^ ,ea:i (C/))(y,n), w/iere is the 

natural extension of the adjoint to the embedding space. 
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